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SECTION  I 
INTRODUCTION 

In  recent  years  there  has  been  a  reawakening  of  interest  in  the  use 
of  ramjet  propulsion  systems  for  modern  strategic  and  tactical  missiles. 
The  advantages  of  ramjet  propulsion  are  the  high  specific  impulse 
(resulting  in  enhanced  range  and/or  speed  capability)  and  a  conceptually 
simple  system  requiring  no  moving  parts  other  than  those  associated  with 
the  fuel  control  system.  The  main  disadvantage  has  been  that  a  booster 
system  is  required  to  provide  the  high  initial  velocity  for  the 
propulsion  system  to  operate. 

The  modern  ramjet  technology  employs  the  Integral  rocket  ramjet 
concept  where  sequential  use  is  made  of  the  ramjet  chamber  as  a  solid 
rocket  combustion  chamber  (for  the  initial  boost  stage  of  flight) 
followed  by  Its  later  use  as  a  ramjet  combustion  chamber  (after  the 
solid  fuel  has  been  burned  and  high  velocity  attained).  Two  conceptions 
for  the  integral  rocket/ramjet  configuration  are  shown  in  Figures  la,  b, 
and  c. 

One  of  the  more  critical  problem  areas  in  the  Integral  rocket/ramjet 
is  the  ramjet  combustion  cycle.  Since  the  ramjet  combustor  must  serve 
the  dual  purpose  as  housing  for  rocket  propellant.  It  cannot  be  ideally 
configured  in  terms  of  fuel  Injection,  flameholders,  and  combustor 
geometry  for  operation  in  the  ramjet  mode.  One  basic  combustor  design  is 
shown  In  Figure  1b.  In  this  type  of  combustor,  fuel  Injection  takes 
place  In  the  air  Inlet  duct  upstream  of  the  dump  station.  Flow 
recirculation  In  the  region  just  downstream  of  the  dump  station  serves  as 
the  primary  flameholding  mechanism.  This  may  be  supplemented  by 
swirling  of  the  inlet  air. 

A  second  combustor  design  is  shown  in  Figure  Id.  In  this 
configuration,  hot  fuel  rich  gas  from  a  solid  fuel  gas  generator  enters 
the  dump  from  a  small  jet  and  mixes  with  air  from  multiple 
circumferentially  located  air  Inlets.  Again,  flow  recirculation 
provides  flame  stabilization. 
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Figure  la.  Integral  Rocket/Ramjet  Missile 
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Figure  Id 
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As  Indicated  in  these  simple  sketches,  the  combustion  flow  is 
multiphase,  turbulent,  and  involves  flow  separation  and  large 
recirculation  regions.  In  addition,  many  ramjet  combustors  of  interest 
contain  fully  3-D  flowfields  because  of  multiple  air  inlets  located  at 
discrete  circumferential  locations  as  illustrated  in  Figures  lc  and  d. 

The  present  work  is  concerned  with  the  computer  modelling  of  the 
ramjet  combustor  flowfields  and  is  part  of  several  concurrent 
theoretical  and  experimental  efforts  undertaken  to  gain  further  insights 
into  the  characteristics  of  these  type  of  flowfields. 
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SECTION  II 
PREVIOUS  WORK 

The  modelling  of  ramjet  combustors  has  been  recently  summarized  by 
Lilley  (Reference  1).  In  his  paper,  lilley  outlines  two  approaches 
which  can  be  distinguished  as  simplified  models  and  Navier-Stokes  type 
solvers.  The  simplified  models  of  the  flow  in  the  past  have  been  very 
popular  (Refrences  2  through  6)  since  they  avoid  the  problem  of  solving 
partial  differential  equations  by  dividing  the  flow  into  regions  such  as 
perfectly  stirred  reactors,  well  stirred  reactors  and  plug  flow  reactors 
with  separate  simple  empirical  models  for  each  region.  These  models  may 
give  useful  qualitative  trends. 

In  recent  years,  however,  there  has  been  a  dramatic  increase  in  the 
popularity  of  methods  based  on  the  direct  numerical  solution  of  the 
Navier-Stokes  equations  for  combusting  flow.  Increases  in  the 
understanding  of  physical  phenomena,  computing  power,  numerical 
algorithm  efficiency,  and  cost  of  experiments  have  all  contributed  to 
this  trend.  The  trend  has  been  the  greatest  in  the  area  of  gas  turbine 
and  furnace  design  and  today  major  engine  manufacturers  use  numerical 
codes  in  their  combustor  design  process  to  predict  performance  trends, 
predict  durability  problem  areas,  and  reduce  testing.  Since  the  codes 
need  to  calculate  the  value  of  all  dependent  variables  in  the  governing 
equations  for  all  points  in  the  flow,  the  distribution  of  a  property  such 
as  density  (for  which  no  suitable  experimental  measurement  method 
presently  exists)  can  be  determined  by  these  models.  Modification  of 
these  codes  for  the  application  to  ramjet  combustor  design  Is  a  feasible 
proposition  (References  7,  8).  The  aim  of  this  report  will  be  to 
outline  present  work  at  AFWAL/PORT  involving  adapting  and  modifying 
developed  numerical  codes  to  the  ramjet  geometry  and  flow  conditions. 
Areas  in  which  future  research  is  needed  will  also  be  outlined. 
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SECTION  III 

THEORETICAL  BACKGROUND 

The  basis  of  any  numerical  model  of  combusting  flow  is  a  set  of  partial 
differential  equations  that  govern  the  flowfield  and  chemistry  of  the 
particular  regime  of  interest.  Formulation  of  these  equations  necessitates 
modelling  of  physical  phenomena.  The  general  structure  of  a  numerical 
code  is  outlined  in  Figure  2. 

There  is  wide  variation  between  codes  regarding  solution  technique 
(Reference  9)  and  the  sophistication  and  extent  to  which  various  physical 
effects  are  included.  The  codes  at  AFVIAL/PORT  include  a  simple  code  for 
modelling  2-0  isothermal  recirculating  flow  and  a  relatively  complex  code 
for  modelling  a  three-dimensional  multiphase  combustor. 

The  codes  in  this  report  all  proceed  in  the  following  manner.  The 
ensemble  averaged  Navier-Stokes  equations  for  the  steady  turbulent  flow 
of  a  viscous  fluid  in  which  turbulent  density  fluctuations  are  ignored 
and  the  transport  equations  for  an  ensemble  averaged  scalar  variable  may 
all  be  written  in  the  following  similar  form  (cylindrical  polar 
coordinates) . 

r c  h (prU4,)  + 1?  (prV*}  *  If  (rr*  !£}  ■  If  (rr*  !£)]  s  s<j>  (1) 

The  term  <p  is  a  general  dependent  variable  and  the  equation  isolates  the 
terms  for  the  convection  of  <J>,  diffusion  (via  turbulent  flux  terms  of  4>) 
and  the  source  (both  creation  and  dissipation)  of  <J>.  The  equations  differ 
not  only  In  their  exchange  coefficient  but  also,  and  primarily.  In 
their  source  term,  S. .  Terms  such  as  pressure  gradients  and  chemical 
reactions  are  treated  In  the  generalized  source  S^. 

To  provide  closure  of  the  equations,  it  is  necessary  to  model  the 
turbulence.  The  method  used  in  this  work  is  the  familiar  k-e  (Reference  10) 
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Figure  2.  Numerical  Code  Structure 
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model  In  which  the  eddy  viscosity  u  can  be  defined  in  terms  of  the  scalar 
quantities,  turbulence  energy  k  and  turbulence  dissipation  rate  e  as  in 
Equation  2. 


For  a  description  of  a  2-D  axisymmetric  isothermal  flow,  the  equations 
may  be  solved  for  <J>  equal  to  the  axial,  radial,  and  swirling  velocities 
(u,v,w),  the  turbulent  kinetic  energy  k  and  turbulent  dissipation  rate  e 
with  the  exchange  coefficients  and  source  terms  outlined  in  Table  B-l . 

Note,  however,  that  the  equations  could  also  be  solved  for  fuel  fraction, 
mixture  fraction,  stagnation  enthalpy,  and  possibly  other  variables  as 
well.  Note  also  that  no  transport  equation  exists  for  pressure  and  that 
the  velocity  terms  must  also  satisfy  the  continuity  equation.  The 
solution  to  this  dilemma  is  to  determine  the  pressure  from  the  Poisson 
equation  (a  combination  of  the  continuity  and  momentum  equations). 

The  technique  used  for  the  solution  of  the  differential  equations  is 
based  on  the  TEACH  program  (References  11,  12).  A  finite  difference 
primitive  variable  (velocities  and  pressure)  procedure  is  used.  An  initial 
guess  is  made  for  pressure  and  the  equations  solved  for  velocity  terms 
on  ASD  CDC  computer.  Corrections  are  then  made  to  velocity  and  pressure 
to  ensure  continuity  and  the  remaining  equations  are  solved.  An  implicit 
line-by-line  relaxation  technique  is  used  in  the  solution  procedure.  A 
staggered  grid  system  is  used  in  which  the  vectors  are  located  midway 
between  the  points  at  which  scalars  are  calculated  (Figure  3).  Walls 
are  simulated  by  the  staircase  method  in  which  the  boundaries  are  drawn 
between  the  points  at  which  the  scalars  are  calculated. 

The  finite  difference  equations  for  each  <j>  are  obtained  by  integrating 
equation  1  over  the  appropriate  control  volume  and  expressing  the  result 
in  terms  of  neighboring  grid  point  values.  The  convection  and  diffusion 
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THREE  GRIDS:  FOR  p,  w  ETC.  -  AT  POSITION  MARKED  (0) 
FOR  u  VELOCITY-  AT  POSITION  MARKED  (-) 
FOR  v  VELOCITY-  AT  POSITION  MARKEO  (t) 

Figure  3.  Staggered  Grid  and  Notation  for  the  Rectangular 
Computational  Mesh 
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terms  are  expressed  as  surface  integrals  and  the  source  term  is 
summarized  resulting  in  an  equation  of  the  form. 


fc**  -  r*  iVe  -  -  r*  IfVw  *  Ep*>  -  r»  UVn  '  [<>V*  '  r*  !&4s  * 

.  .  (3) 

[Sp<p  +  Sj]  X  Vol 

To  enhance  convergence,  only  terms  which  are  negative  are  permitted  in  Sp. 


Boundary  conditions  are  problem  dependent;  however,  the  following 
general  conditions  apply.  Inlet  values  are  directly  specified.  At  the 
outlet,  axial  velocities  are  calculated  to  ensure  mass  continuity,  radial 
velocities  are  set  to  zero  and  zero  normal  gradient  is  specified  for 
remaining  terms.  Near  wall  tangential  velocities  are  connected  to  their 
zero  wall  values  by  way  of  tangential  shear  stress  wall  functions.  Near 
wall  values  for  e  are  fixed  using  length  scales  near  the  wall  and  the 
current  value  of  k. 
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SECTION  IV 

PRESENT  COOES  AT  AFWAL/PORT 

There  are  at  present  five  numerical  codes  at  AFWAL/PORT  and  these  are 
briefly  outlined  below. 

1 .  TEACH-X 

This  is  a  basic  version  of  the  Imperial  College  Teach  Code  suitable 
for  isothermal,  2-dimensional  recirculating  turbulent  flows.  The  code 
solves  u  and  v  velocity  components,  pressure,  turbulence  kinetic  energy 
and  turbulence  dissipation  in  isothermal,  incompressible  flow.  The 
program  input  is  designed  for  dump  combustors. 

2.  TEACH- T 

This  code  is  similar  to  TEACH-X  but  includes  two  extra  equations; 
one  for  the  transport  of  a  second  species,  and  one  for  the  transport  of 
specific  enthalpy.  The  mixture  fraction  and  temperature  (determined 
from  the  enthalpy)  in  conjunction  with  the  equation  of  state  determine 
the  density  fluctuations  in  the  flowfield.  Note  that  the  fluid  kinetic 
energy  terms  have  been  neglected  in  the  enthalpy  equation.  The  program 
includes  combustion  only  to  the  extent  that  a  single  product  is  produced 
by  Instantaneous  ignition  and  infinitely  rapid  reaction  kinetics.  This 
first  approximation  is  only  suitable  for  diffusion  limited  combustion 
where  there  Is  no  initial  mixing  of  fuel  and  air.  The  program  Inputs 
are  designed  for  a  dump  type  combustor  equipped  with  a  concentric  annulus 
burner  at  the  dump  plane. 

3.  STARPIC 

The  title  of  this  code  Is  an  acronym  for  swirling  turbulent 
axisymmetric  recirculating  flows  in  practical  isothermal  geometries.  The 
code  was  prepared  by  D.  G.  Li 1 1 ey  and  0.  L.  Rhode  of  Oklahoma  State 
University  under  NASA  grant  NAG-3-74  (Reference  7).  This  code  is 
similar  to  the  Teach-X  code  previously  described  except  that  an  equation 
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for  a  swirl  velocity  component  is  included.  The  code  was  written 
essentially  to  model  the  flow  in  turbine  combustors  and  the  subroutine 
that  specifies  the  initial  conditions  is  somewhat  more  flexible  than  the 
TEACH  codes.  The  inlet  conditions  may  include  a  sloping  upstream 
boundary  which  is  modelled  within  the  calculation  regime  by  a  series  of 
boundary  steps. 

Provision  is  also  made  to  read  the  initial  conditions  from  a  previous 
solution  tape.  This  procedure  is  handy  but  does  not  fulfill  the  role  of 
a  true  restart  capability  since  it  is  necessary  ta  run  the  solution  to 
completion  to  obtain  the  solution  tape;  and  if  the  job  is  aborted  through 
some  calculation  error  or  job  time  limit,  the  solution  derived  by  the 
program  at  that  stage  is  lost. 

The  program  also  includes  routines  for  calculation  of  non-dimensional 
values,  stream  functions,  and  plot  files.  These  are  useful  calculations 
but  would  be  better  included  in  a  post  processor  since  they  are 
extremely  wasteful  of  memory  space. 

The  program  contains  variable  underelaxation.  At  present,  this  is 
in  a  somewhat  primitive  form  in  that  the  underelaxation  factors  increase 
linearly  between  fixed  values  according  to  some  preset  conditions 
determined  by  the  programmer.  Winterfeldt  has  indicated  that  considerably 
improved  convergence  can  be  achieved  by  relating  the  underelaxation 
factor  to  actual  values  of  the  local  variables  and  geometry  (Reference  13). 
This  Is  In  accordance  with  the  general  experience  of  high  underelaxation 
factors  leading  to  divergence  of  the  solution  and  low  underelaxation  factors 
leading  to  excessively  long  convergence  times.  Nevertheless,  It  should 
be  stated  that  the  use  of  even  the  present  method  of  varying  the  underelax¬ 
ation  factor  which  places  considerable  reliance  on  the  experience  of  the 
programmer  is  superior  to  the  retention  of  fixed  underelaxation  factors. 

4.  STARRC 

This  title  is  acronynmous  for  Swirling  Turbulent  Axlsymmetrlc 
Recirculating  and  Reacting  Compressible  Code.  This  code  was  written  by  the 
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present  author  using  the  aforementioned  STARPIC  as  a  base  code  to  which 
several  features  were  added.  The  diffusion  equation  from  TEACH-T  was 
added  together  with  the  enthalpy  equation  except  that  the  enthalpy  could 
now  be  either  the  specific  static  or  stagnation  enthalpy  as  determined 
by  a  logic  switch.  The  inclusion  of  the  equation  of  state  for  calculation 
of  density  together  with  the  stagnation  enthalpy  equation  allows 
isentropic  compressibility  effects  to  be  included.  Note,  however,  the 
solution  of  the  pressure  equation  has  not  yet  been  altered  to  include 
variations  in  density  so  that  the  solution  tends  to  be  explicit  rather 
than  implicit  with -resultant  convergence  problems.  The  boundary  conditions 
allow  the  inclusion  of  a  sloping  upstream  boundary  and  a  downstream  exit 
nozzle  both  of  which  are  modelled  by  a  series  of  finite  steps.  An 
expanding  or  an  expanding-contracting  grid  can  automatically  be  generated 
in  the  x  direction.  The  first  approximation  combustion  of  the  TEACH-T 
code  has  been  included  and  provision  made  for  the  planned  addition  of 
reaction  kinetics.  A  post  processor  for  producing  various  2-dimensional 
and  3-dimensional  plots  has  been  written  by  Schwartzkopf  for  this  code. 

The  format  for  inputting  initial  conditions  is  somewhat  more  general 
than  that  of  the  STARPIC  code  on  ASD  CDC  computer. 

The  entire  code  has  been  written  as  a  series  of  correction  decks  for 
Lilley's  STARPIC  Code  on  ASD  CDC  computer. 

5.  THE  3-D  CODE  PERFORMANCE  MODEL  PROGRAM 

This  program  was  prepared  for  NASA  under  contract  NAS3-22542  by 
Garrett  Turbine  Engine  Company  (Reference  14)  and  is  the  extension  of  an 
earlier  code  produced  for  the  US  Army  and  described  in  Report  No  USARTL- 
TR-55C.  The  3-D  program  Is  general  and  capable  of  predicting  recirculating 
flows  In  3-dimensional  geometries.  At  present,  the  code  Inputs  and 
boundary  conditions  are  specifically  oriented  towards  gas  turbine 
combustor  geometries.  Reacting  or  non-reacting,  swirling  or  non-swirling, 
diffusion  and/or  premixed  flames,  and  gaseous  and/or  liquid  fuel  combustion 
can  be  handled  by  the  program.  The  code  solves  for  all  the  quantities 
of  STARRC  in  a  3-dimensional  flowfield  plus  the  mass  fractions  of 
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unburned  fuel,  oxygen,  carbon  monoxide,  C  H  _  (the  intermediate  Hydro- 

x  yt 

carbon  in  a  four  step  reaction)  H^,  C02  and  H^O,  three  radiation  flux 
vectors,  soot  and  NO  emissions  and  the  fuel  spray  trajectory,  droplet 
size  distribution,  and  evaporation  rate. 


The  program  includes  the  following  physical  models: 

a.  Turbulence  -  the  two  equation  (k-e)  turbulence  model 

b.  Chemistry  -  A  four  step  chemical  reaction  scheme  as  outlined  below 


CxHy 

CxHy-2  +  I  °2 
CO  +  >s0, 


CxHy-2  +  H2 


xCO  +  ^  H2 


CO. 


H2  +  »s02 


H20 


c.  Chemical  Reaction  Rate 

The  reaction  rates  are  governed  by  the  minimum  of  either  the  time 
averaged  Arrhenius  model  of  the  reaction  kinetics  or  the  component  mixing 
rate  as  determined  by  a  turbulent  eddy  break  up  model . 


d.  Soot  Emissions 

The  program  contains  a  quasiglobal  model  that  requires  the  solution 
of  transport  equations  for  the  concentration  of  nuclei  and  soot.  Two 
particle  sizes  for  soot  are  assumed,  a  small  size  resulting  from 
nucleatlon  and  a  large  size  resulting  from  fuel  droplet  pyrolosls  and 
char  formation. 


e.  Radiation 

A  six  flux  model  Is  used.  The  absorption  coefficient  Is  calculated 
locally  as  a  function  of  concentration  of  soot,  carbon  dioxide,  and 
water  vapor. 
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f.  NOx  emissions 

NOx  emissions  are  calculated  using  the  CREK  kinetic  model. 

A  typical  solution  procedure  using  the  Garrett  Code  involves  three 
steps  executed  automatically  in  the  following  manner: 

(1)  Solution  of  all  required  variables  except  soot  radiation 
and  NOx  till  cumulative  mass  residual  of  approximately  5%  is  reached. 

(2)  Inclusion  of  soot  and  radiation  in  the  solution  procedure 
and  continue  till  cumulative  mass  residual  of  1%  is  reached. 

(3)  Inclusion  of  NOx  equations  and  continue  solution  until  the 
desired  convergence  level  for  the  final  solution  (approx.  0.5$)  is 
reached . 

The  present  application  of  thes’e  five  codes  is  directed  primarily 
towards  the  use  of  STARRC  and  the  NASA-Garrett  3-0  Code.  STARRC  is 
suitable  for  axially  symmetric  cold  flows  and  is  configured  for  operation 
on  the  AFWAL  CDC  7600  computer  for  geometries  of  Figure  1c.  The  NASA- 
Garrett  Code  is  suitable  for  3-dimensional  geometries  and  is  presently 
being  configured  for  operation  on  both  the  AFWAL  computer  and  the  Cray 
computer  at  Klrtland  AFB  for  geometries  of  Figure  Id.  Both  of  these 
codes  are  receiving  further  development  and  will  be  further  detailed 
in  the  following  sections. 

At  the  time  of  writing  this  report,  a  sixth  code  Is  being  prepared  to 
model  isothermal  flow  in  the  three  dimensional  geometries  of  Figure  Id. 
Since  this  code  will  not  contain  the  dependent  variables  associated  with 
combustion  processes.  It  will  be  far  more  conservative  in  memory 
requirements  and  consequently  will  allow  greater  resolution  than  the 
NASA-Garrett  code.  This  code  will  be  supplementary  to  the  NASA-Garrett 
code  and  will  principally  be  used  for  the  modelling  of  water  tunnel  and 
cold  gas  flows  in  the  geometries  of  Figure  Id. 
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SECTION  V 

THE  COMPUTER  CODE  STARRC 


1 .  INTRODUCTION 

As  mentioned  previously,  the  STARRC  code  draws  heavily  on  Lilley's 
STARPIC  Code  (Reference  7)  and  has  in  fact  been  written  as  a  series  of 
UPDATE  correction  decks- to  the  STARPIC  Code.  These  correction  decks  will 
be  outlined  and  the  following  sections  should  be  read  in  conjunction  with 
Lilley's  report. 

The  present  program  library  version  of  STARRC  incorporates  all  these 
correction  sets  in  a  resequenced  format  with  the  exception  of  the  common 
block  correction  sets.  The  correction  sets  are  discussed  individually 
simply  to  outline  their  individual  reasons  for  inclusion  and  are  not 
separately  identified  in  the  resequenced  format. 

A  brief  description  of  the  operation  of  STARRC  and  the  STARRC  post 
processor  (intended  as  a  supplementary  users  manual)  and  a  summary  of 
present  experience  follows  the  correction  deck  descriptions. 

2.  CORRECTION  DECK  DESCRIPTIONS 

a.  NEWCOMMON  9624,  N0N0N0IM4824 

These  two  correction  decks  modify  the  common  blocks  and  the  definition 
of  IT,  JT  (Indices  of  maximum  dimensions  of  dependent  variables)  to  allow 
the  array  structures  to  accommodate  the  required  grid  size.  NEWCOMMON  9624 
allows  a  96  x  24  grid  and  simple  editing  of  values  96  and  24  will  allow 
any  desired  grid  size.  N0N0NDIM4824  modifies  to  a  48  x  24  grid  and 
removes  the  nondlmenslonal  and  stream  function  calculation  routines. 

This  enables  a  48  x  24  grid  problem  to  be  run  within  the  day  shift 
memory  restriction  of  170,000g  words. 
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b.  Ident  DREWGRID 

This  correction  set  applies  the  boundary  conditions  for  one  case  of 
the  DREWRY  (Reference  15)  experiment  discussed  (Case  1). 

c.  Ident  CALCFR 

This  correction  deck  solves  the  finite  difference  equation  for  the 
transport  of  mixture  fractions. 

d.  Ident  CALCH 

This  correction  deck  solves  the  finite  difference  equation  for  the 
transport  of  enthalpy. 

e.  Ident  COMBUST 

This  correction  deck  modified  the  main  controlling  subroutine  and 
property  calculation  subroutine  to  include  the  calculations  of  density, 
mixture  fraction,  and  enthalpy.  Density  may  be  fixed  for  the  case  of 
fixed  density  isothermal  flow  or  calculated  from  the  equation  of  state 
using  the  local  values  of  absolute  pressure,  mixture  fraction  and 
temperature.  If  stagnation  enthalpy  is  calculated  (INMACH  *  .TRUE,  and 
1NCALCH  *  .TRUE.),  compressibility  effects  will  be  included  since  local 
temperature,  and  hence  local  density,  will  be  determined  from  local 
enthal py. 

f.  Ident  COMCOM 

This  correction  set  modified  the  common  blocks  to  Include  the  enthalpy 
and  fuel  mixture  fraction  terms. 

g.  Ident  COMPRO 

This  correction  set  modified  the  subroutine  PROMOO  to  include 
chapters  for  enthalpy  and  fuel  mixture  fraction. 
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h.  Ident  CORECT 

This  set  provides  minor  corrections  to  Lilley's  original  code 

i.  Ident  ENERGY 

This  correction  set  alters  the  way  in  which  the  zero  axial  gradient 
condition  is  applied  at  the  exit  plane.  The  enthalpy  values  at  the 
last  axial  grid  points  are  gradually  altered  from  initial  values  to  ensure 
conservation  of  energy.  This  modification  was  included  in  an  attempt  to 
improve  the  convergence  for  nonadiabatic  combustor  walls  in  future 
calculations . 

j.  Ident  INLET 

This  ident  provides  code  to  automatically  generate  a  stepped  grid 
for  a  sloping  inlet  wall.  The  Input  parameter  is  inlet  axial  length, 
ALINLT.  The  code  generates  steps  about  a  curve  defining  the  expansion 
of  the  wall  radius  from  inlet  radius  to  combustor  radius.  The  general 
curve  describing  the  wall  radius  is,  at  present,  a  linear  expansion; 
however,  other  curves  may  easily  be  substituted.  This  ident  made  obsolete 
a  previous  correction  deck  ENTRAN  which  generated  a  dump  inlet  only. 

k.  Ident  EXIT 

This  ident  provides  code  to  automatically  generate  a  stepped  grid  for 
simulation  of  an  exit  nozzle.  The  code  also  supplies  the  boundary 
conditions  necessary  for  an  exit  nozzle.  There  are  two  input  parameters, 
ALTRAN  and  JEXIT.  ALTRAN  is  the  axial  distance  from  the  combustor  inlet 
to  the  point  at  which  the  transition  to  the  radial  reduction  of  the  nozzle 
occurs.  JEXIT  is  the  J  index  of  the  grid  point  located  immediately  inside 
the  exit  radius.  The  code  generates  steps  about  a  curve  describing  the 
radius  of  the  nozzle  as  a  function  of  axial  location.  Two  general  curves 
have  been  programmed  -  a  linear  contraction  and  a  convex  elliptical 
contraction. 
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l.  Idents  -  ROGRID,  ROGFLO,  RUNDAT 

These  correction  decks  provided  the  grid  dimensions,  flow  data,  and 
run  data  for  a  typical  test  case. 

m.  I  dent  NEWSWIRL 

This  ident  was  written  to  include  inlet  swirl  in  the  form  of  a  free 
vortex  profile  used  in  the  experiments  of  Buckley  et  al  (Reference  16). 
Note  the  redefinition  of  the  parameter  NSBR  for  defining  type  of  swirl. 

A  new  parameter  is  RHUB.  This  determines  the  radius  within  which  the 
inlet  swirl -is  considered  to  have  solid  body  rotation  irrespective  of 
the  nature  of  the  remainder  of  the  swirl  profile.  Note  the  definition 
of  inlet  solid  body; swirl  velocity  in  terms  of  input  swirl  number  as 
supplied  by  Li  1  ley ’{Equation  4).  This  contrasts  with  the  definition 
by  Buckley  (Equation  5)  though  the  methods  of  the  authors  for  calculation 
of  swirl  number  for  a  given  flow  are  identical 

(4) 

(5) 

SWNB  *  Input  Swirl  Number 


u  _  SWNB  *  *  R 

W  "  T *SWNB  U  Rstep 


W  =  SWNB  *  U  * 


Rstep 


n.  Ident  REPAIR 

This  ident  contains  miscellaneous  corrections  to  the  code.  These 
have  been  written  In  response  to  problems  discovered  in  the  operation  of 
the  code.  These  problems  may  have  resulted  from  inclusion  of  other 
correction  sets  or  simply  from  the  way  the  initial  conditions  were 
specified. 

o.  Ident  TITLE 

This  ident  added  conment  cards. 
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p.  Ident  OPGRID,  CHGRID 

These  correction  decks  enable  the  optional  selection  of  a  grid  spacing 
that  geometrically  increases  in  the  axial  direction  away  from  the  inlet 
plane,  or  alternately,  a  grid  spacing  that  geometrically  increases  from 
both  the  inlet  plane  and  the  nozzle  exit  plane  towards  the  center  of  the 
combustor . 

q.  Obsolete  Idents  may  appear  in  the  list  of  idents  in  the  program 
library.  These  obsolete  idents  include  MODIFY  and  XXX. 

3.  OPERATIONAL  DETAILS  OF  THE  STARRC  CODE 

The  STARRC  code  contains  16  subroutines.  The  flow  chart  for  the 
program  is  shown  in  Figure  4.  The  functions  of  the  individual  subroutines 
can  be  summarized  as  follows: 

a.  MAIN 

Controls  and  monitors  the  entire  sequence  of  calculations:  initial¬ 
ization,  properties  and  initial  output;  the  iteration  loop  with  calls  to 
update  main  variables,  other  mixture  properties  and  intermediate  output; 
and,  after  termination  of  the  iteration  loop,  final  output,  an  increment 
in  inlet  degree  of  swirl  and  a  return  to  the  beginning  again. 

b.  INIT 

Sets  values  to  the  numerous  geometric  quantities  concerned  with  grid 
structures,  and  initalizes  most  variables  to  zero  or  other  reference 
value. 


c.  PROPS 

Updates  the  fluid  properties  via  sequential  calculation  of  mass 
fractions  of  fuel  oxidizer  and  products,  temperature,  density  and 
turbulent  viscosity.  Uses  underrelaxation  for  the  last  two  of  these,  the 
k-e  turbulence  model  and  appeals  to  PROMOD  for  any  other  modifications. 
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MAIN  CONTROLLING 
SUBROUTINE 
■CONTRO' 


INITIAL  |  MAIN 

CONDITIONS  LOOP 


Figure  4.  STARRC  Flowchart 
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d.  PRINT 

Prints  out  an  entire  variable  field  according  to  a  standard  format. 

e.  CALCU  and  CALCV 

Calculates  coupl ing  coefficients  of  finite  difference  equation  for 
axial  velocity  u*  and  radial  velocity  v*.  calls  PROMOD  for  boundary 
modifications  and  LISOLV  for  entire  field  of  variables  to  be  updated  to 
get  u*  and  v*  fields. 

f.  CALCP 

Calculates  coupling  efficients  of  finite  difference  equation  for 
pressure  correction  p';  calls  PROMOD  for  boundary  modifications  and 
LISOLV  to  obtain  p'  field.  The  subroutine  closes  with  p*,  u*  and  v* 
being  'corrected'  with  p',  u',  and  v'. 

g.  CALCW,  CALCTE,  CALCED,  CALCFR,  CALCH 

Calculates  coupling  coefficients  of  appropriate  finite  difference 
equation,  calls  appropriate  part  of  PROMOO  and  then  LISOLV  for  complete 
update  of  the  variable  in  question. 

h.  PROMOD 

Modifies  the  values  of  the  finite  difference  equation  coefficients, 
or  the  variables,  near  walls  or  other  boundaries  where  particular 
conditions  apply.  The  subroutine  is  divided  into  chapters,  each  handling 
a  particular  variable  and  being  called  from  a  CALC  subroutine,  and  each 
chapter  considers  all  the  boundaries  around  the  solution  domain. 

i.  LISOLV 

Updates  entire  field  of  a  particular  variable,  by  applying  TOMA 
(tridiagonal  matrix  algorithm)  to  all  the  lines  in  the  r-dlrection 
sequentially  from  left  to  right  of  the  integration  domain. 
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From  this  description  it  can  be  seen  that  in  general  the  user  is 
interested  in  making  modifications  to  only  3  subroutines;  CONTRO,  PROPS 
and  PROMOD  with  the  majority  of  corrections  being  applied  to  subroutine 
CONTRO.  In  the  present  format  the  code  does  not  read  any  input.  All 
input  data  (other  than  restart  data)  is  contained  as  Fortran  statements 
and  recompilation  is  necessary  for  each  execution.  The  normal  operating 
procedure  is  therefore  to  write  the  input  code  as  a  correction  set  to 
the  basic  code.  This  premise  is  used  in  the  -operating  instructions  set 
out  below. 

The  following  instructions  are  referenced  to  subroutine  name  and  line 
number  (STARRC,  CY  =15)  and  are  intended  as  an  operational  checklist  for 
the  novice  to  use  in  conjunction  with  Lilley's  report  (Reference  7). 

1.  Set  Common  Block  Sizes  and  values  for  IT  and  JT  using  NEWCOMMONnn 
or  NONONDIftnm . 

2.  CONTRO. 88  CONTRO. 89 

Data  block  entries  for  swirl  numbers  and/or  vane  blade  angles  for  the 
solution  of  problems  with  swirl.  The  program  always  assumes  that  swirl 
is  present  so  that  to  solve  for  only  the  zero  swirl  case  it  is  necessary 
to  specify  a  single  solution  and  set  the  first  parameters  of  the  swirl 
data  block  equal  to  zero. 

3.  CONTRO. 96 

This  statement  deletes  error  messages  resulting  from  exponent 
underflow  and  sets  the  variable  for  which  the  error  was  detected  to  zero. 

4.  CONTRO. 98  C0NTR0.118 

Logical  switches.  The  majority  of  these  are  self  explanatory, 
however  the  following  rules  may  be  helpful. 

a.  With  the  present  parameters,  I  FINE  *  .FALSE,  implies  a  coarse 
grid  of  48  points  in  the  axial  direction.  Of  these  48  points,  45  are 
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generated  to  cover  the  specified  geometry  and  3  (set  by  NEXRAC)  are  set 
for  constant  diameter  geometry  downstream  of  the  exit.  The  axial  locations 
for  these  last  3  points  are  obtained  by  reflecting  the  previous  3  points 
about  the  exit  plane.  If  I  FINE  *  .TRUE.,  a  fine  grid  is  generated  with 
96  points  in  the  axial  direction.  Of  these  96  points,  91  are  located 
within  the  specified  geometry,  and  5  (set  by  NEXRAF)  are  set  for  constant 
diameter  geometry  downstream  of  the  exit  and  their  axial  location  is 
determined  by  reflection  of  the  previous  5  points  about  the  exit  plane. 

b.  INCONT  =  .TRUE,  implies  that  the  grid  spacing  in  the  axial 
direction  expands  to  the  center  of  the  combustor  and  then  contracts  at 
the  same  rate  to  the  exit  plane  of  the  combustor.  If  INCONT  is  set  to 
.FALSE.,  the  grid  spacing  will  continue  to  expand  in  the  axial  direction 
to  the  exit  plane  of  the  combustor. 

c.  INMACH  *  .TRUE,  implies  calculation  of  stagnation  enthalpy  rather 
than  static  enthalpy. 

d.  INREAC  *  .TRUE,  implies  chemical  reaction  occurs. 

e.  IROTEM  *  .FALSE,  implies  fixed  density  isothermal  flow. 

IROTEM  *  .TRUE,  implies  calculation  of  density  from  equation  of  state. 

If  INMACH,  IROTEM  and  INCALH  are  all  set  .TRUE,  compressibility  effects 
are  included  in  the  flow  calculations. 

5.  CONTROL. 122  -  see  instruction  52. 

6.  CONTROL. 124  -  NSTLN  *  n.  The  value  n  sets  the  number  of  streamlines 

calculated.  The  default  value  is  11. 

7.  CONTRO.125  -  NPLTLN  *  n.  The  value  n  sets  the  number  of  streamlines 

plotted  on  the  line  printer  plot.  The  default  value  is  6. 

8.  C0NTR0.126  -  MAXLN  *  n.  The  value  n  sets  the  maximum  number  of 
streamlines  that  may  be  plotted.  The  default  value  is  10. 
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9.  CONTRO . 1 28  -  JPRINT  =  NITER  +  n.  The  value  n  determines  the  number 
iterations  between  printing  of  field  variable  values.  The  default  value 
is  1100. 

10.  CONTRO. 129  -  IPRINT  =  NITER  +  n.  The  value  n  determines  the  number 
of  iterations  between  printing  residual  sums  and  monitor  values  of  field 
variables.  The  default  value  is  1. 

11.  CONTRO. 130  -  LFS  =  n.  The  value  n  sets  the  index  of  the  vane  blade 
angle  array  or  the  index  of  the  swirl  number  array  for  the  first  swirl 
calculation.  The  default  value  is  1. 

12.  CONTRO. 131  -  LFS MAX  =  n.  The  value  n  sets  the  index  of  the  vane 

blade  angle  array  or  the  index  of  the  swirl  number  array  for  the  last 

swirl  calculation.  The  default  value  is  1. 

13.  CONTRO. 133  -  MAXIT  «  NITER  +  n.  The  value  n  is  the  maximum  number 

of  iterations  to  be  run  for  each  swirl  case  if  the  solution  is  not  stopped 
by  either  a  convergence  or  a  divergence  criterion.  The  default  value  is 
1000. 

14.  CONTRO. 132  -  NSBR  »  n.  The  value  n  determines  the  type  of  swirl 

profile.  Three  values  for  n  are  presently  recognized. 

a.  n  *  1  Solid  Body  Rotation  from  Swirl  Generator 

b.  n  =>  2  Flat  Swirl  Profile  from  Swirl  Vanes 

c.  n  =  3  Free  Vortex  Swirl  Profile 

The  present  default  value  is  1. 

15.  CONTRO. 135  CONTRO. 136  -  IT  3  m  JT  3  n.  See  instruction  1. 

16.  CONTRO. 138  CONTRO. 145  -  NSWP0  3  n0  -  The  value  n0  sets  the  number 
of  application  of  the  line  iteration  for  the  solution  of  the  0  equation. 
Present  devault  values  are 


a.  NSWPU  3  4 
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b.  NSW  Pi*  =  5 

c.  Remaining  NSWP0  =  3 

The  general  recommendation  is  that  NSWPU  should  be  greater  than  other 
velocities  and  scalars  with  the  exception  of  P.  NSWPP  should  be 
approximately  twice  NSWP0.  If  the  problem  of  divergence  arises,  the 
generally  recommended  solution  is  the  increase  of  NSWP0  and  the  decrease 
of  URF0.  The  present  author  has  found  little  effect  results  from  the 
increase  of  NSWP0  with  the  exception  of  NSWPP.  In  problems  that  exhibit 
divergence,  this  may  be  increased  to,  say,  NSWPP  =  15  with  some  small 
effect  on  increase  in  stability. 

17.  C0NT0.192  -  ISTEP  *  n.  The  value  n  sets  the  axial  location  of  the 
first  node  inside  the  inlet  plane.  The  default  value  is  2. 

18.  C0NT0.193  -  JSTEP  *  n.  The  value  n  sets  the  radial  location  of  the 
first  node  interior  to  the  inlet  wall.  The  default  value  is  14. 

19.  C0NTR0.194  -  JEXIT  =  n.  The  value  n  sets  the  radial  location  of 
first  node  interior  to  the  exit  nozzle  wall  at,  or  downstream  of,  the 
exit  plane.  The  default  value  is  14. 

20.  CONTRO.199  -  NJ  *  n.  The  value  n  sets  the  total  number  of  nodes  in 
the  radial  direction  for  which  values  are  assigned  or  calculated.  The 
default  value  Is  24. 

21.  C0NT0.205  -  R LARGE  *  X.  The  value  X  is  the  radius  of  the  combustor 
The  default  value  is  0.048768. 

22.  C0NTR0.198  -  INDCOS  *  n.  The  value  n  defines  the  coordinate  system 
n  =  1  plane  flows,  n  *  2  axisymmetric  flows. 

23.  C0NTR0 . 206  -  ALTOT  *  X.  The  value  X  defines  the  length  of  the 
combustor  from  the  inlet  plane  to  the  exit  plane.  The  default  value  Is 
0.4064  m. 
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24.  C0NTR0.207  -  ALTRAN  =  X.  The  value  n  defines  the  length  of  the 
combustor  from  the  inlet  plane  to  the  transition  point  at  which  the  area 
reduction  of  the  nozzle  commences.  The  default  value  is  0.381  m. 

25.  C0NTR0.208  -  ALINLT  =  X.  The  value  X  defines  the  length  of  the 
combustor  from  the  inlet  plane  to  the  transition  point  where  constant  or 
maximum  cross  sectional  area  occurs.  The  default  value  is  0. 

26.  C0NTR0.21 1  C0NTR0.212  -  NI  =  m  NEXRAC  =  n.  See  note  4a,  m  and  n 
are  settings  for  coarse  mesh. 

27.  C0NTR0.214  -  EPSX  =  X.  The  value  X  sets  the  geometric  ratio  for 
axial  grid  spacing  increase/decrease  for  a  coarse  mesh.  The  default 
value  is  1.11. 

28.  C0NTR0.250  C0NTR0.251  -  NI  *  m  NEXRAF  *  n.  See  note  4a,  m  and  n 
are  settings  for  fine  mesh. 

29.  C0NTR0.253  -  EPSX  *  X.  The  value  X  sets  the  geometric  ratio  for 
axial  grid  spacing  increase/decrease  for  a  fine  mesh.  The  default  value 
is  1.102. 

30.  C0NTR0.289  -  Y(j)  *  X.  The  value  X  is  the  radial  location  of  the 
jth  node  in  the  radial  direction. 

31.  C0NTR0.331  -  RINLT  «  f(x(I)).  The  function  f(x(I ) )  Is  the  definition 
of  the  Inlet  expansion  profile.  Present  default  is  linear  expansion. 
Alternate  maximum  radial  nodes  are  located  either  side  of  defined  profile. 

32.  C0NTR0.366  -  RTRAN  *  f(X( I ) )  -  The  function  f(x(I))  is  the  definition 
of  the  nozzle  contraction  profile.  Present  default  Is  linear  contraction. 

33.  C0NTR0.382  C0NTR0.391  -  INCAL0  =  .TRUE.  Defining  the  value  INCAL0 

to  be  TRUE/FALSE  allows/inhlbits  the  solution  of  the  equation  for 
variable  <t>. 
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34.  C0NTR0.397  CONTRO. 400  -  C  =  X.  The  value  of  constants  In  the  k-e 

3 

turbulence  model  are  defined  here. 

35.  C0NTR0.404  C0NTR0.407  -  PR0  -  X.  Turbulent  Prandtl /Schmidt  numbers 

are  defined  here.  Default  values  are  PRED  *  1.22,  PRTE  =  1.0,  PRH  =  0.9, 

PRFR  =0.9. 

36.  C0NTR0.409  -  UIN  =  X.  The  value  X  is  the  average  inlet  velocity. 

This  value  is  used  in  determining  the  initial  conditions. 

37.  CONTROL. 411  -  TURBIN  =  X.  The  value  of  n  is  the  turbulence  intensity 
level  at  the  input.  The  default  value  is  2%. 

38.  C0NTR0.413  -  ALAMDA  =  X.  The  value  of  internally  defined  length 
scales  is  X  times  the  appropriate  length.  The  initial  boundary  value 
for  turbulence  dissipation  is  determined  from  a  length  scale 
ALAMDA*C0MB  US  TOR  DIAMETER. 

39.  C0NTR0.415  -  VISCOS  =  X.  Laminar  Viscosity  is  defined  by  X.  The 
present  default  value  is  17.11  x  10"®  kg/ms. 

40.  C0NTR0.417  -  PRANDT  =  X.  Laminar  Prandtl  number  is  equal  to  X.  The 
present  default  is  0.7. 

41.  CONTRO.  418  -  TEMP  *  n.  Gas  inlet  temperature  is  defined  by  X. 

The  present  default  value  is  248. 5°K. 

42.  CONTRO. 420  -  HFU  ■  X.  The  fuel  heat  of  reaction  is  defined  by  X. 

The  present  default  value  is  0.,  representing  the  situation  of  no  combustion. 

43.  CONTRO. 422  -  CPR  »  X.  The  quantity  X  represents  the  average 
specific  heat  at  constant  pressure  for  the  gas  mixture.  The  default 
value  is  1004. 
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44.  C0NTR0.424  -  WPR  =  X  WOX  =  Y  WFU  *  Z.  The  quantities  X,  Y,  Z, 
represent  the  average  molecular  weights  of  combustion  products,  oxidizer, 
and  gaseous  fuel  respectively.  The  default  values  are  given  for  air  as 
the  oxidizer  and  for  fuel  simulated  by  argon. 

45.  C0NTR0.428  -  XI  =  X.  The  value  X  is  the  stochiometric  oxidizer/fuel 
mass  ratio. 

46.  C0NTR0.432  -  FUEL  =  X.  The  value  X  is  the  fuel  mass  fraction 
entering  the  combustor. 

47.  C0NTR0.436  C0NTR0.438  -  PO  =  X  IPREF  =  m  JPREF  =  n.  The  value  and 

location  of  the  reference  pressure  is  defined  here.  The  values  of  m,  n 
should  be  specified  within  the  calculation  domain.  The  default  location 
is  (2,2),  near  the  intersection  of  symmetry  axis  and  the  inlet  plane. 

48.  C0NTR0.440  C0NTR0.442  -  IMON  =  m  JMON  =  n  SORMAX  =  X.  The  quantities 

m,  n  defined  the  node  for  which  the  monitor  values  of  the  three  velocities, 
pressure  and  turbulence  dissipation  rate  are  printed  at  intervals 
determined  by  IPRINT.  The  value  defined  by  SORMAX  is  the  maximum  accepted 
value  for  the  sum  of  residual  sources  for  convergence  to  occur. 

49.  C0NTR0.477  -  U(2,n)  *  X,  TE(l,n)  =  Y,  ED(l,n)  =  Z.  Axial  velocity, 
turbulence  kinetic  energy  and  turbulence  dissipation  rate  at  the  inlet 
plane  for  all  radial  grid  points  are  defined  here. 

50.  C0NTR0.507  -  RHUB  *  X.  The  radius  of  the  hub  of  the  swirler  is 
defined  here.  The  hub  radius  is  the  region  within  which  the  swirl  at  the 
Inlet  plane  is  assumed  to  be  solid  body  rotation. 

51.  C0NTR0.534  -  T(1,J)  =  X.  Inlet  plane  temperature  profile  defined 
here  equal  to  X. 

52.  C0NTR0.570  -  P(1,J)  a  X.  The  Inlet  plane  pressure  profile  defined 
here  equal  to  X.  Default  value  is  flat  profile  equal  to  reference 
pressure. 
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53.  C0NTR0.693  -  Write  (11)  LIST.  This  statement  writes  to  the  restart 
tape  a  series  of  values  to  be  used  for  normalization  of  dependent 
variables  in  the  output  plots.  When  rereading  the  restart  tape,  these 
values  are  read  into  the  array  DIMLES. 

54.  C0NTR0.663,  664  C0NTR0.740  C0NTR0.754,  761  -  URF0  =  n.  This 

section  sets  the  underrelaxation  factors.  The  default  values  should 
prove  satisfactory  for  constant  density  flow.  Some  adjustment  may  be 
necessary  for  varying  density  flow.  Decreasing  the  underrelaxation 
factor  improves  stability  at  the  expense  of  convergence  rate.  For 
varying  density  flows,  it  is  best  to  start  with  the  fixed  density 
solution  as  the  initial  condition.  Another  tip  is  to  use  a  LO-HI-LO 
sequence  for  underrelaxation  factors  to  bring  the  solution  to  quick  stable 
convergence.  Further  work  is  needed  in  this  area. 

4.  OPERATION  OF  STARRC  POST  PROCESSOR 

The  present  STARRC  Post  Processor  is  named  TCHPLT.  The  code  will 
generate  3-0  perspective  views  or  2-0  plots  of  families  of  curves  of  the 
general  variable  Z  versus  the  X  and  Y  coordinates  using  the  OISSPLA 
package. 

To  operate  TCHPLT,  it  is  necessary  to  first  generate  (IWRITE  »  .TRUE.) 
a  restart  file  in  STARRC.  This  file  has  the  local  file  name  TAPE11  and 
the  normal  procedure  is  to  catalog  this  file  at  the  end  of  a  STARRC  run. 
This  file  is  then  reattached  as  TAPE9  for  TCHPLT. 

TCHPLT  is  run  in  batch  mode  with  the  job  control  cards  of  Appendix  A 
and  a  data  deck  constructed  in  the  following  manner.  All  data  fields 
are  10  characters  in  width. 

a.  3-Dimensional  Plot 

Card  1:  (A10,  110,  6A10)  3DPL0T,  DATAFILE,  PLOT  TITLES  -  The  string 

3DPL0T  is  left  justified.  The  DATAFILE  Is  a  right  justified  Integer 
defining  the  location  of  the  file  on  the  tape  of  the  dependent  variable 
(Table  B-2)  to  be  plotted  on  the  Z  axis.  If  datafile  is  a  negative 
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number  -k,  all  quantities  on  file  k  will  be  normalized  by  the  values 
written  on  the  file  by  STARRC  (See  5.3.53).  PLOT  TITLE  is  a  user-defined 
plot  label  and  must  be  terminated  with  a  %. 

Card  2:  (8110)  XDIM,  FIRSTX,  LASTX  -  These  three  quantities  are 

right  justified  integers  defining  the  dimension  of  the  X  array  and  the 
indices  of  the  first  and  last  X  locations  to  be  plotted  respectively. 

Card  3:  (8110)  YDIM,  FIRSTY,  LASTY  -  These  three  quantities  are 

right  justified  integers  defining  the  dimension  of  the  Y  array  and  the 
indices  of  the  first  and  last  Y  location  to  be  plotted  respectively. 

Card  4:  (2E10.0,  6A10)  XAXIS,  XNORM,  XLABEL? 

Card  5:  (2E10.0,  6A10)  YAXIS,  YNORM,  YLABEL? 

Card  6:  (2E10.0,  6A10)  ZAXIS,  ZNORM,  ZLABEL?  -  The  values  XAXIS,  YAXIS 

ZAXIS,  define  the  proportional  lengths  of  the  respective  AXES.  The  values 
XNORM,  YNORM,  ZNORM  may  be  any  positive  real  number  and  are  normalization 
factors  for  the  appropriate  axes.  The  ALPHA  strings  XLABEL,  YLABEL,  ZLABEL 
provide  the  labels  for  the  appropriate  axes  and  must  be  terminated  by  %. 

Card  7:  BLANK  CARD 

Card  8:  (2E10.0,  6A10)  PHI,  THETA  -  PHI  and  THETA  define  the  view 

angle.  PHI  is  the  angle  positive  anticlockwise  from  XAXIS,  THETA  is  the 
angle  positive  when  the  observer  Is  above  the  XY-plane. 

b.  Two-Dimensional  Plots 

Card  1:  (A10,  110,  6A10)  2DLABEL,  DATAFILE,  PLOT  TITLE?  -  The  alpha 

string  2DLABEL  determines  the  type  of  plots  according  to  the  following 
table. 
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LABEL  PLOT 

2DYVSX  PLOTS  of  Y  vs  X  for  constant  Z 

2DZVSX  PLOTS  of  Z  vs  X  for  constant  Y 

2DZVSY  PLOTS  of  Z  vs  Y  for  constant  X 

The  value  of  Z  is  the  value  of  the  variable  defined  by  the  DATAFILE. 
Comments  for  DATAFILE  and  PLOT  TITLE  for  3-dimensional  plots  apply. 

Cards  2,  3,  4,  5,  6:  These  cards  are  identical  to  those  used  in  3-D 
plots  with  the  exception  that  the  values  XAXIS,  YAXIS,  ZAXIS  now  define 
the  actual  length  of  the  plot  AXES  in  inches. 

Card  7:  (8E10.0)  XMIN,  DELX,  YMIN,  OELY  -  This  card  sets  the  plot 

scaling.  A  blank  card  will  invoke  automatic  scaling.  To  predetermine 
the  scales,  set  XMIN  and  YMIN  to  minimum  value  on  X  and  Y  axes  respectively 
and  DELX  and  DELY  to  the  required  scale,  in  units  per  inch  of  plot  for  X 
and  Y  axes  respectively. 

Card  8:  (8110)  N  -  The  quantity  N  is  a  right  justified  integer 

defining  the  number  of  curves  in  the  family  to  be  plotted. 

Card  9:  (8110)  a-j ,  c^.  a3»  an  "  ^rd  8  contains  N  values,  one  for 

each  requested  curve  in  the  family  plotted.  The  term  defines  the  value 
of  the  parameter  held  constant  for  the  ith  curve. 

Card  9  to  9+N:  (E10.0,  2110)  S(I),  JJS(I),  IIS(I)  -  Cards  9  to  9+N 

contain  parameters  for  the  N  requested  curves.  The  term  S(I)  sets  the 
value  of  the  parameter  held  constant  for  the  ith  curve.  The  term  JJS(I) 
sets  the  symbol  to  be  used  in  plotting.  JJS(I)  may  have  a  value  between 
0  and  15.  The  term  I IS ( I )  determines  whether  symbols  are  plotted 
according  to  the  following  table 

ISYM  =  -1  point  symbols  only 

ISYM  *  0  connecting  lines  only 

ISYM  *  +1  point  symbols  and  connecting  lines 
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c.  Termination  of  Data  Deck 

The  data  deck  must  be  terminated  with  a  card  containing  the  left 
justified  alpha  string  DONE. 

5.  OPERATIONAL  EXPERIENCE  WITH  STARRC  CODE 

This  section  outlines  a  number  of  problems  encountered  with  the 
STARRC  Code  and  possible  solutions. 

a.  Iteration  Control 

The  iteration  process  is  monitored  by  comparison  of  the  absolute 
values  of  the  residual  sources  of  mass  and  preselected  dependent  variables 
in  the  flowfield  with  a  present  value,  SORMAX,  representing  the  maximum 
source.  Iteration  is  terminated  when  the  largest  residual  source  is  less 
than  SORMAX,  a  divergence  criteria  is  invoked,  or  a  preset  number  of 
iterations  is  reached.  In  general,  this  procedure  has  worked  well  and 
the  approach  is  particularly  suited  for  production  runs  from  restart 
tapes  where  an  approximate  idea  exists  of  the  number  of  iterations  needed 
for  solution. 

When  a  new  case  is  being  modelled  for  the  first  time,  the  procedure 
is  not  so  successful.  With  this  type  of  problem,  interactive  programming 
is  preferred  to  allow  problem  debugging,  monitoring  of  convergence,  and 
alteration  of  quantities  such  as  underelaxation  factors  as  required. 
Unfortunately,  core  memory  restrictions  prevent  this  on  the  present 
computer  system. 

A  proposed  solution  Is  to  remove  the  iteration  number  limit  and 
Instead  to  relay  on  job  time  limit  for  maximum  limit  to  the  number  of 
iterations.  After  each  iteration,  the  execution  time  (from  subroutine 
second  (t)  on  CDC  system)  is  compared  with  the  job  time  limit  and  at  a 
preset  difference  the  restart  file  is  generated  and  the  execution 
terminated.  This  will  prevent  the  loss  of  the  restart  file  that  presently 
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occurs  if  job  time  limit  is  reached  before  convergence  or  iteration 
limit  and  will  allow  the  batch  submission  of  a  number  of  short  time  jobs 
to  hasten  debugging  and  monitor  convergence. 

b.  Convergence  Control 

The  problems  of  convergence  are  essentially  those  of  excessive 
slowness  to  converge,  divergence,  or  convergence  with  lack  of  conservation 
of  variables  that  require  conservation.  The  Fortran  variables  which 
influence  the  iteration  behavior  are  the  number  of  update  sweeps  NSWP0 
and  the  under-relaxation  factors  URF0.  As  stated  earlier,  if  divergence 
occurs,  the  remedy  generally  lies  in  increasing  the  former  (especially 
for  pressure)  and  decreasing  the  latter.  For  excessively  slow  convergence, 
the  reverse  applies.  Note,  however,  that  increasing  the  number  of  grid 
points  also  tends  to  increase  the  number  of  iterations  required  for 
convergence.  In  general,  the  present  values  of  NSWPJ9  and  the  routine 
adapted  from  STARPIC  for  varying  URF0  have  been  found  adequate  for 
constant  density  flows. 

The  solution  is  not  so  clear  for  varying  density  flows.  For  many 
problems  convergence  is  still  obtained  simply  by  increasing  NSWPP  to  15, 
reducing  URFP  to  0.5  and  URFRHO  to  0.2,  rather  starting  the  varying  density 
solution  with  initial  conditions  derived  from  a  wholly  or  partly  converged 
constant  density  solution.  The  reasons  that  this  procedure  is  not  always 
successful  are  inherent  in  the  method  in  which  the  density  variations 
have  been  introduced.  Consider  the  effect  of  the  various  logic  switches. 

If  IROTEM  is  set  .TRUE.,  local  density  is  calculated  from  the  equation 
of  state  and  rather  than  assigned  a  fixed  initial  value.  Variations  in 
local  absolute  pressure  will  alter  the  local  density.  If  INCALFR  is  set 
.TRUE.,  a  second  species  is  transported  and  local  density  variations  may 
also  occur  due  to  variations  in  concentration  of  a  gas  of  different 
molecular  weight.  If  INCALH  is  set  .TRUE.,  local  variations  in  temperature 
and,  through  the  equation  of  state,  density  may  occur  since  temperature 
is  Tiow  calculated  from  local  enthalpy  rather  than  assigned  a  fixed  value. 

If  INMACH  is  set  .FALSE,  and  INCALH  is  set  .TRUE.,  the  static  enthalpy  is 
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calculated  and  temperature  variations  will  occur  only  if  heat  is  released 
due  to  chemical  reaction  ( INREAC  = .TRUE . ) .  If  however,  INMACH  is  set  .TRUE, 
and  INCALH  is  set  .TRUE.,  stagnation  enthalpy  is  calculated  and  variations 
in  local  kinetic  energy  will  be  reflected  in  variations  in  temperature, 
density,  and  pressure. 

Under  all  of  the  above  conditions,  variations  in  density  arise  only 
through  the  equation  of  state  in  the  PROP  subroutine.  The  equation  for 
calculation  of  pressure  has  not  been  altered.  Consequently,  the  solution 
proceeds  in  the  following  way. 

(1)  With  an  assumed  pressure  p*  and  density  p*  initial  values 
of  axial  and  radial  velocities  u*  and  v*  are  calculated. 

(2)  Using  the  above  values  and  the  equations  of  continuity  of 
mass  and  momentum,  corrected  values  u,  v,  and  p  are  derived  u  =  u*  +  u', 
v  =  v *  +  v 1 ,  p  =  p*  +  p'. 

(3)  Using  the  new  values  for  p  and  the  latest  value  of  temperature 
a  new  value  for  p  is  defined  from  the  equation  of  state. 

As  a  result  of  step  (3)  above,  local  continuity  will  no  longer  occur 
and  the  solution  for  pressure  and  density  tends  to  be  explicit  rather 
than  implicit  with  possible  divergence.  If  the  value  of  pressure  used 
in  the  calculation  of  the  equation  of  state  is  large  compared  to  the 
pressure  correction  terms,  density  correction  terms  are  small  and 
convergence  results.  Otherwise,  the  solution  will  tend  to  diverge. 
Increasing  the  grid  spacing  changed  a  diverged  solution  to  a  convergent 
solution  in  one  test  case.  This  procedure  is  not  always  acceptable  and 
the  long  term  solution  is  to  alter  the  algorithm  for  calculation  of 
pressure  correction  terms  to  include  the  effect  of  variations  in  density. 

If  the  solution  procedure  stagnates  or  converges  without  conservations 
of  quantities  that  should  be  conserved,  present  experience  suggests  that 
this  nearly  always  results  from  an  inconsistent  set  of  boundary  conditions. 
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Krishnamurthy  (Reference  17)  reports  that  a  coarse  grid  mash  can  also 
lend  to  an  absence  of  mass  conservation  for  the  diffusion  of  a  second 
species . 

c.  Accuracy  of  Solution 

The  accuracy  with  which  the  differential  equations  are  being  solved 
is  dependent  on  a  number  of  factors,  one  of  which  is  grid  size.  If 
increases  in  mesh  density  produce  no  significant  changes  in  the  dependent 
variables  at  locations  in  the  flow  of  interest  to  the  modeller,  the  grid 
configuration  may  be  considered  as  adequate.  The  grid  densities  used  in 
the  following  case  studies  may  be  used  as  guides. 

It  is  interesting  to  note  here  the  problem  that  occurs  in  the 
specification  of  the  step  boundary  condition  for  both  sloping  inlets  and 
nozzles.  The  problem  occurs  if  a  step  has  several  points  in  the  horizontal 
direction  and  is  manifested  as  an  oscillating  rather  than  steady  change 
in  the  static  pressure  along  the  wall.  This  problem  is  illustrated  in 
Figure  5  which  shows  the  configuration  of  a  stepped  nozzle  and  indicates 
the  experimentally  computed  static  pressures  (relative  to  the  flow 
reference  pressure).  It  is  doubtful  if  these  variations  have  any 
significant  effect  on  the  flowfield  and  it  is  assumed  that  the  problem 
is  numerical.  Some  further  Investigation  of  the  boundary  conditions  for 
the  solution  of  the  pressure  correction  equation  might  prove  helpful. 

The  simulation  of  the  nozzle  boundary  by  a  nonorthogonal  grid  (thus 
eliminating  the  present  staircase  procedure)  should  remove  this  problem. 

To  a  large  extent,  assessment  of  the  accuracy  of  the  modelling 
procedure  for  othe"*  than  simple  pipe  flows  suffers  from  a  dearth  of  data 
from  well  defined  experimental  programs  that  could  be  used  as  a  baseline 
for  comparison  purposes.  Similarly,  predictions  of  flows  for  combustor 
geometries  is  often  hampered  by  lack  of  knowledge  of  boundary  values. 
Several  experimental  programs  are  planned  to  correct  these  deficiencies. 

Three  case  studies  follow  that  illustrate  the  use  of  the  program  at 
the  present  stage  of  development. 
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6.  CASE  STUDIES 

Case  1  Model  Verification  for  Ax i symmetric  Dump  Combustor 
Geometry 

Though  a  number  of  research  programs  are  attempting  to  remedy  the 
situation,  at  present  one  of  the  most  comprehensive  investigations  of 
the  turbulent  flow  in  an  axisymmetric  dump  combustor  geometry  is  that  of 
Drewry  (Reference  15)  who  has  investigated  cold  flow  conditions  employing 
surface  flow  visual ization,  pressure  probing  and  on-line  gas  sampling 
with  a  quadrupole  mass  spectrometer .  This  case  study  reports  the  modelling 
of  one  of  Drewry's  configurations  illustrated  in  Figure  6.  An  array  of 
eight  circumferential  fuel  injection  ports  located  63.5  mm  upstream  of 
the  dump  plane  was  connected  to  a  separate  gas  (fuel)  feed  system. 
Measurements  included  wall  static  pressure,  static  and  total  pressure 
traverses  and  concentration  measurements  made  when  argon  was  injected 
through  the  fuel  ports.  It  is  worth  noting  that  the  concentration 
measurements  were  performed  in  a  horizontal  plane  (two  fuel  injectors 
are  in  this  plane)  and  no  check  was  made  on  the  axis  of  symmetry  of  the 
measured  fuel  distribution. 

Drewry's  first  traverse  station  for  both  velocity  and  concentration 
was  made  at  2.54  mmm  from  the  dump  plane.  This  traverse  was  used  to 
set  the  initial  conditions  of  velocity  and  fuel  mass  fraction.  A 
uniform  inlet  static  pressure  profile  and  2%  turbulence  intensity  were 
assumed.  The  combustor  was  modelled  by  a  48  x  24  variably  spaced  grid 
for  the  variable  density  solution  and  a  96  x  24  variably  spaced  grid  for 
the  fixed  density  solution.  Rapid  convergence  occurred  for  the  fixed 
density  solution.  The  variable  density  solution  was  more  sensitive  to 
initial  conditions  and  chosen  underrelaxation  factors  and  convergence 
was  not  achieved  for  the  96  x  24  grid. 

Measured  and  predicted  values  for  the  wall  static  pressure,  velocity 
distribution  and  fuel  mass  fraction  are  shown  in  Figures  7,  8,  and  9. 

Note  that  excellent  agreement  occurs  between  predicted  and  measured  values 
of  wall  pressure  and  velocity.  The  difference  between  the  model  and  the 
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experiment  for  the  velocity  in  the  region  of  the  wall  can  be  attributed 
to  the  fact  that  Drewry  was  using  pitot  probes  in  a  highly  turbulent 
region  of  the  flow  and  that  the  model  is  using  stationary  axisymmetric 
flow  while  the  experiments  observed  that  this  was  not  the  true  case.  Note 
that  the  coarser  grid  variable  density  model  gave  slightly  better 
predictions  of  velocity  and  fuel  mass  fraction  and  slightly  worse 
predictions  of  wall  pressure  than  the  finer  grid  fixed  density  model. 

The  model  needs  to  be  refined  to  enable  use  of  variable  density  without 
the  special  care  presently  required  to  obtain  convergence. 

The  fuel  distribution  shows  good  qualitative  trends  but  there  is 
disagreement  between  the  actual  and  predicted  values.  The  normalized 
profiles  would  show  considerably  better  agreement.  The  difference  can  be 
attributed  to  the  fact  that  the  fuel  distribution  was  not  axisymmetric 
prior  to  the  dump  plane  and  that  the  input  profile  used  had  higher  average 
than  the  actual  average  value  (since  it  was  in  the  plane  of  2  fuel 
injectors  as  previously  indicated).  The  present  plan  is  to  further  extend 
the  verification  of  the  computer  model  by  comparison  of  other  dump 
configurations  of  Drewry  (Reference  15)  and  of  Boray  (Reference  18). 

Case  Study  2  -  Effect  on  Fuel  Distribution  of  Addition  of  Swirl 

Buckley  et  al  (Reference  16)  have  shown  that  significant  reduction  in 
combustion  chamber  length  and  total  pressure  loss  can  be  achieved  by  the 
introduction  of  swirl  to  the  flow  upstream  of  the  dump  plane.  This  case 
study  investigates  the  distribution  of  the  total  fuel  throughout  the 
combustor  chamber  in  the  flow  conditions  of  case  1  when  swirl  is  introduced. 
Three  swirl  profiles  (Forced  vortex  or  solid  body  rotation.  Flat  swirl 
or  fixed  angle,  and  Free  vortex)  were  investigated. 

Results  are  presented  for  nominal  swirl  number  0.4  for  the  extremes  of 

Forced  Vortex  (W^r)  and  Free  Vortex  (VK-)  swirl  -  the  Flat  swirl  results 

r 

lie  between  these  extremes. 
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For  all  profiles,  the  inlet  dump  plane  contained  a  region  of  forced 
vortex  (solid  body  rotation)  from  the  centerline  to  .1875  Rstep.  This 
represented  the  hub  used  in  the  swirlers  produced  by  Buckley  (Reference  16). 
Calculations  were  made  for  48  x  24  grid  variable  density  flow.  Figure  10 
shows  the  fuel  mass  fraction  distribution  throughout  the  combustor  for 
no  swirl,  forced  vortex  and  free  vortex  swirl. 

Note  that  the  free  vortex  profile  produces  an  evenly  distributed  fuel 
mass  fraction  more  rapidly  than  the  forced  vortex.  Qualitatively,  these 
results  show  agreement  with  the  experiments  of .Buckley  who  found  that 
the  swirl  profiles  could  be  ranked  Free  vortex,  Constant  angle,  or  Forced 
vortex  in  order  of  decreasing  combustion  efficiency.  Care  must  be  taken 
not  to  extrapolate  this  comparison  too  far  since  even  the  forced  swirl 
case  gave  a  better  combustion  efficiency  (Reference  16)  than  no  swirl 
while  this  was  not  the  case  for  the  computer  predictions  of  the  fuel 
mass  fraction  distribution.  The  reason  for  the  variation  with  swirl 
profile  lies  principally  in  the  effect  of  the  swirl  pattern  on  the 
central  and  peripheral  recirculation  zones.  It  should  be  noted  that 
swirl  profile  is  only  one  of  several  factors  that  may  effect  these 
recirculation  zones.  Other  factors  are  reported  by  Rhode  and  Li 1 1 ey 
(Reference  1 9) . 

Several  authors  (References  19-22)  have  attempted  to  model  swirl  in 
combustor  configurations  and  have  encountered  problems  with  the  prediction 
of  turbulence.  It  is  prudent  therefore  to  recommend  further  experiments 
Involving  swirl  with  the  measurement  of  turbulence  and  species  in  both 
hot  and  cold  flows.  These  experiments  are  under  consideration  at 
AFWAL/PORT. 

Case  3  -  Modifications  to  Outer  Wall 

This  case  study  relates  to  the  alteration  of  combustor  outer  wall 
geometry  in  an  effort  to  reduce  the  total  pressure  loss  through  the 
combustor.  The  initial  configuration  is  shown  in  Figure  11.  The 
combustor  has  a  shorter  L/D  than  that  of  Drewry  and  contains  an  elliptical 
rather  than  a  straight  converging  nozzle. 
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The  proposed  configuration  change  suggested  by  Buckley  was  to  utilize 
a  sloping  wall  rather  than  a  dump,  thus  eliminating  the  peripheral 
recirculation  zone  and  hopefully  reducing  total  pressure  loss.  The  flame 
stabilization  was  to  be  provided  by  a  swirl  induced  central  recirculation 
zone.  The  combustor  was  modelled  for  the  flow  of  cold  air  (no  fuel)  using 
a  48  x  24  variably  spaced  grid  and  the  variable  density  solution.  Free 
vortex  swirl  (nominal  swirl  number  0.4,  0.6)  was  introduced  for  the 
proposed  configuration  change  of  45°  sloping  inlet  and  constant  angle 
inlet  to  the  start  of  the  nozzle.  The  flowfield  streamlines  for  the 
latter  geometry  and  for  a  dump  combustor  for  the  case  of  no  swirl  and 
nominal  swirl  number  0.4  are  depicted  in  Figure  11.  Examination  of 
Figure  11  reveals  that  the  peripheral  recirculation  zone  is  removed  by 
the  introduction  of  the  sloping  wall  at  the  expense  of  a  reduction  in 
size  of  the  central  recirculation  zone.  Originally  it  was  hoped  that  the 
central  recirculation  zone  would  increase  in  size  with  the  alteration  in 
geometry.  This  reduction  in  the  central  recirculation  zone  could  lead 
to  flame  stability  problems  at  reasonable  swirl  numbers  and  consequently 
no  combustion  experiments  were  performed  with  the  altered  geometry. 

An  interesting  proposed  follow-on  case  study  would  be  to  run  the 
code  with  the  inlet  fuel  distribution  of  Drewry  and  observe  the  effect 
on  the  fuel  distribution  in  cold  flow  of  these  configuration  changes. 
Experimental  measurement  could  be  included  with  those  mentioned  in  case 
study  2  to  see  if  the  central  recirculation  zone  is  actually  decreased 
and  what  effect  this  has  on  the  fuel  distribution.  The  results  of  these 
numerical  and  experimental  studies  would  be  a  better  guide  than  the  results 
presently  available  to  the  wisdom  of  proceeding  with  combustion 
experiments. 
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SECTION  VI 

THE  NASA-GARRET  CODE 


1.  INTRODUCTION 

The  plan  for  the  adaptation  of  this  code  to  ramjet  configurations 
proposes  a  number  of  separate  tasks  as  outlined  below. 

a.  Adapt  code  for  operation  on  ASD  CDC  computer. 

b.  Adapt  code  for  operation  on  AFWL  Kirtland  CRAY  computer. 

c.  Adapt  code  for  side  inlet  combustor  with  no  dome  inlet  flow.  Model 
this  configuration  for  water  flow  and  compare  results  with  water  tunnel 
measurements. 

d.  Add  a  dome  inlet  flow  to  the  model  of  case  3  to  simulate  the 
presence  of  the  gas  generator  jet  in  the  combustor.  Model  the  config¬ 
uration  for  water  flow  and  compare  results  with  that  of  the  water  tunnel 
experiments  and  the  results  predicted  by  a  3-0  isothermal  flow  code 
capable  of  finer  grid  resolution  than  the  NASA-Garrett  code  (Section  IV). 

e.  Model  the  flow  of  cold  air  and  compare  results  with  those 
obtained  in  gas  sampling  experiments  and  Laser  Doppler  Anemometer 
measurements  and  the  results  predicted  by  a  3-0  isothermal  flow  code 
capable  of  finer  grid  resolution  than  the  NASA-Garrett  code  (Section  IV). 

f.  Model  the  combustion  of  Hydrocarbon  fuel  utilizing  geometries 
presently  under  experimental  evaluation.  Compare  results  with  experiments. 

g.  Perform  experiments  with  the  code  Involving  variations  of  present 
experimental  parameters.  This  may  Involve  variations  In  geometry,  inlet 
flow  conditions,  and  the  nature  of  the  fuel. 

2.  PROGRESS  TO  DATE 

Tasks  a,  b,  c  have  been  proceeding  simultaneously  and  are  reported 
collectively.  The  principal  reason  for  adapting  the  code  to  the  CRAY 
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computer  is  to  make  use  of  its  superior  memory  storage  and  speed 
capabilities.  Approximately  an  order  of  magnitude  decrease  in  the 
runtime  can  be  achieved  without  major  changes  in  the  code  logic  by 
running  the  code  on  the  CRAY  vector  processor  rather  than  the  ASD  CDC 
machine.  The  present  system  for  submission  of  jobs  to  the  CRAY  from  ASD 
is  not  an  acceptable  procedure  for  debugging  of  programs  because  of 
limited  access  availability  and  slow  response.  Consequently,  the  best 
procedure  is  to  debug  codes  on  the  ASD  machine  and  to  limit  CRAY 
operations  to  production  runs.  Job  submission  to  the  CRAY  may  include 
the  program  and/or  data  as  input  or  may  retrieve  these  from  the  mass 
storage  system  at  Kirtland. 

Operation  of  the  code  on  the  ASD  CDC  computer  is  straightforward. 

The  code  must  be  compiled  in  Fortran  5  because  of  the  nature  of  the 
format  statements.  The  principal  limitation  to  effective  operation  is 
restriction  on  memory  size.  The  supplied  three  dimensional  code  requires 
an  execution  file  length  of  367,200g  words  for  a  non-optimized  array 
structure  suitable  for  a  15  x  15  x  7  grid.  This  requirement  is  too  close 
to  the  imposed  third  shift  limit  of  377,000g  words  to  allow  loading  of 
the  debug  routines.  Consequently,  a  general  update  correction  deck 
named  "MYCOMMONBLOCKUPDATER"  has  been  constructed  to  allow  optimum  array 
construction  for  the  desired  grid  size.  "MYCOMMONBLOCKUPDATER"  is  edited 
according  to  instructions  supplied  in  Appendix  B  and  the  resultant  file, 
for  examDle  "COMBUG",  is  applied  as  an  UPDATE  correction  deck  to  the 
"THREEDPROGRAMLIBRARY" .  "COMBUG"  will  allow  the  construction  of  a 
15x7x7  grid  and  when  loaded  with  debug  options  requires  313,400g 
words  of  memory.  Further  Increases  In  grid  size  may  be  achieved  by 
removing  the  calculations  of  Nox  emissions  from  the  program.  Modifications 
to  the  code  allowing  the  simulation  of  the  side  entry  combustor  are 
included  in  the  update  correction  deck  "SIDEINUPDATER" .  The  procedure 
has  been  to  simulate  the  side  Inlet  by  altering  the  boundary  conditions 
for  the  locations  designated  as  radial  injection  points  in  Reference  14 
to  provide  the  injection  velocity  with  a  radial  and  an  axial  component. 

This  procedure  is  relatively  unsophisticated  and  may  need  to  be  Improved. 
The  remainder  of  the  correction  deck  Is  to  allow  for  the  conditions  of 
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zero  fuel,  zero  dome  inlet  flow,  and  constant  density  flow  as  required 
for  the  comparison  with  the  water  tunnel  experiments.  It  should  also  be 
noted  that  the  following  corrections  should  be  made  to  the  present  card 
input  deck: 

a.  Card  20,  fuel  properties  must  contain  non-zero  values  even  if 
no  fuel  is  input. 

b.  Card  28,  last  three  fields  should  read  TCYLW,  TINLW.  TUP. 

c.  Card  26,  description  of  dome  inlet.  JSW1 ,  JSW2  must  contain 
values  greater  than  1  even  if  no  flow  through  this  inlet. 

d.  Card  27,  temperature  at  the  dome  inlet  must  not  be  set  at  zero 
even  if  no  flow  through  this  inlet. 

Further  modification  to  this  deck  are  being  made  to  enhance  "user 
friendliness"  of  this  code.  It  will  also  be  necessary  to  supply  fuel 
injection  at  the  side  inlets  so  that  the  combustion  experiments  of 
task  6  may  be  simulated. 
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SECTION  VII 

CONCLUSION  AND  FUTURE  WORK 

The  current  design  practice  for  simple  gas  turbine  and  industrial 
combustor  geometries  employs  computer  codes  in  preference  to  experi¬ 
mentation  for  predicting  parametric  trends.  Future  development  of  these 
codes  will  obviously  lead  to  wider  application  for  more  complex  flows. 

The  present  work  has  illustrated  that  the  codes  presently  under 
development  for  gas  turbine  combustors  can  also  be  applied  to  ramjet 
combustor  configurations.  Qualitative  trends  can  be  predicted  by 
present  codes  and  the  determination  of  the  quantitative  accuracy  of  the 
codes  is  somewhat  limited  by  the  lack  of  well  defined  experimental  data. 

Future  work  in  the  ramjet  modelling  area  may  proceed  along  a  number 
of  paths.  The  following  tasks  can  be  clearly  identified. 

1.  Further  refinement  of  the  present  2-D  model  and  the  expansion  of 
the  model  to  include  reaction  kinetics  in  order  to  obtain  a  more  realistic 
model  of  the  ramjet  dump  combustion  process.  Further  numerical 
simulations  and  experimental  measurements  are  needed  to  determine  the 
limits  of  applicability  of  the  present  modelling  procedure.  Possible 
experiments  are  suggested  In  Section  V.6. 

2.  A  continuation  of  the  task  of  adapting  the  3-D  codes  to  the 
Ramjet  configurations  of  Figure  Id. 

3.  The  adaptation  of  the  codes  to  a  more  user  friendly  format 
including  such  features  as  grid  generation  and  initial  condition 
preprocessor.  Iteration  diagnostics  and  graphic  displays.  The  generation 
of  a  data  base  of  solutions  on  restart  tapes  to  form  initial  conditions 
for  future  similar  problems  would  be  a  useful  procedure. 
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4.  The  modification  of  the  codes  to  include  such  operational 
improvements  as  mesh  imbedding,  adaptive  grids,  models  of  sloping 
boundaries  and  dynamic  variation  of  underrelaxation  factors. 

5.  Improvement  in  models  of  physical  phenomena.  Consideration 
should  be  given  to  higher  order  closure  models  presently  being  included 
in  some  codes.  The  use  of  solid  fuels  may  necessitate  the  derivation 
of  combustion  reaction  models  peculiar  to  the  ramjet  area. 
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APPENDIX  A 

Typical  TCHPLT  Job  Control  Deck 

JOB,  T25,  1050,  CM  155000.  PXXXXXX. 

ATTACH,  TCHPLT,  STARRCPLOTTER,  ID  =  HARCH. 

ATTACH,  TAPE9,  STARRCWRITEFIIE,  ID  =  XXXXX. 

REQUEST,  TAPE99 ,  *Q. 

ATTACH,  DISSPLA8,  ID  =  A780283,  SN  =  ASD. 

LIBRARY,  DISSPLA. 

FTN,  I  =  TCHPLT. 

MODE,  1  . 

LGO. 

ATTACH,  UNP  1038,  ID  =  LIBRARY,  SN  =  ASD. 

UNPL038. 

ROUTE,  TAPE99,  DC  *  PU,  TID  =  XX,  ST  =  CSA,  FC  *  NG. 

*EOR 

DATA 

*END  OF  JOB 

Note  that  the  FC  parameter  controls  the  plot  paper  according  to  the 
following  format: 

NG  Narrow  Grid 

MG  Wide  Grid 
NP  Narrow  Plain 

WP  Wide  Plain 
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APPENOIX  B 

Editing  of  MYCOMMONBLOCKUPDATER 


An  update  deck  that  will  modify  the  library  deck  "THREEDPROGRAMLIBRARY" 
to  enable  calculations  for  a  grid  L,  M,  N  can  be  produced  by  editing  the 
general  correction  deck  "MYCOMMONBLOCKUPDATER"  in  the  following  manner: 


REPLACE  NXYZ 

7*XYZ 
NX ,NY ,NZ 
NX, NY 
NXY 
M2XYZ 
XYMZ 
NX 
NY 
NZ 

1 XYZP1 
2XYZP1 
3XYZP1 
4XYZP1 
5XYZP1 
6XYZP1 


BY  L*M*N 

7*L*M*N 

L.M.N 

L,M 

L*M 

(L-2)*(M-2)*(N-2) 

(L-2)*(M-2) 

L 

M 

N 

1 *L*M*N+1 
2*L*M*N+1 
3*L*M*N+1 
4*L*M*N+1 
5*L*M*N+1 
6*L*M*N+1 


TABLE  B-l 


THE  FORM  OF  THE  COMPONENTS  OF  THE  LINEARIZED  SOURCE  TERM*, 
THE  CELL  VOLUME  INTEGRAL  /  S^dV  3  S*  <l>p  +  S  j 

OF  EQ.  (1)  FOR  2-D  AXISYMMETRIC  ISOTHERMAL  FLOW 


4» 

r* 

s^/v 

V/v 

1 

0 

0 

0 

u 

U 

0 

<.u  3P 
s  "  3x 

V 

M 

-2  p2 

i/> 

i 

w 

U 

0 

_  £VW  W  3_  {  ) 

r  ^7  3r 

k 

u/ak 

-CyCDp2k/y 

G 

e 

u/oe 

-C2pe/k 

C*!  CyGpk/y 

In  this  table  certain  quantities  are  defined  as  follows: 

-u  9  ,  Buu  1  3  .  3v  x 
S  ’  37  (M  37)+F  3?  (nj  W 


*In  this  table,  V  stands  for  the  cell  control  volume  and  u  «  u 

j 
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G  =  u[2{(^)2  +  (^)2  +  (l)2}  +  (|i  +  |j)2 

+  {rl7^>2  +  ^ 

TABLE  B-2 


FILE  NUMBERS  FOR  TCHPLT  ROUTINE 
LFS  IS  THE  SWIRL  LOOP  INDEX 


File  Number 

Dependent  Variable 

(LFS-1  )*1 3  + 

1 

Axial  Velocity  U 

(LFS-1  )*1 3  + 

2 

Radial  Velocity  V 

(LFS-1)  *13  + 

3 

Swirl  Velocity  W 

(LFS-1)  *13  + 

4 

Static  Pressure  P 

(LFS-1)  *13  + 

5 

Turbulence  Kinetic  Energy  TE 

(LFS-1  )*1 3  + 

6 

Turbulence  Dissipation  ED 

(LFS-1  )*1  3  + 

7 

Effective  Viscosity  VIS 

(LFS-1  )*13  + 

8 

Stream  Function  STFN 

(LFS-1)  *13  + 

9 

Mixture  (total  fuel)  Fraction  FR 

(LFS-1)  *13  + 

10 

Oxidizer  Mass  Fraction  OX 

(LFS-1  )*13  + 

11 

Fuel  Mass  Fraction  FU 

(LFS-1  )*13  + 

12 

Enthal py  H 

(LFS-1)  *13  + 

13 

Density  p 
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